Load libraries

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(tinytable))
suppressPackageStartupMessages(library(rairtable))
suppressPackageStartupMessages(library(lme4))
suppressPackageStartupMessages(library(broom))

Load data

Sample metadata

sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
    mutate(Extraction=factor(Extraction, levels=c("ZYMO","DREX","EHEX"))) %>%
    rename(dataset=Dataset)
extraction_data <- airtable("tblBcTZcRG1E9wsGO", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
  read_airtable(., fields = c("EX DNA ng","Datasets_flat"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
  rename(id=1,extract=2,dataset=3) %>%
  filter(dataset %in% sample_metadata$dataset) %>%
  select(dataset,extract)
library_data <- airtable("tblo6AuYpxbbGw9gh", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
  read_airtable(., fields = c("Required PCR cycles","Datasets_flat"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
  rename(id=1,pcr=2,dataset=3) %>%
  filter(dataset %in% sample_metadata$dataset) %>%
  select(dataset,pcr)
indexing_data <- airtable("tblhfsiR4NI9XJQG0", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
  read_airtable(., fields = c("Adaptors (nM)","Library (nM)","Datasets_flat"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
  rename(id=1,adaptors=2,library=3,dataset=4) %>%
  filter(dataset %in% sample_metadata$dataset) %>%
  mutate(ratio=library/(adaptors+library)) %>%
  select(dataset,adaptors,library,ratio)
preprocessing_data <- airtable("tblJfLRU2FIVz37Y1", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
  read_airtable(., fields = c("bases_pre_fastp","bases_post_fastp","host_bases","metagenomic_bases","singlem_fraction","C","EHI_plaintext"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
  rename(id=1,nonpareil=C,dataset=EHI_plaintext,microbial_fraction=singlem_fraction) %>%
  filter(dataset %in% sample_metadata$dataset) %>%
  select(dataset,bases_pre_fastp,bases_post_fastp,bases_post_fastp,metagenomic_bases,microbial_fraction,nonpareil)
assembly_data <- airtable("tblG6ZIvkYN844I97", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
  read_airtable(., fields = c("EHI_number_api","assembly_length","N50","L50","num_contigs","num_bins","metagenomic_bases"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
  rename(id=1,assembly_mapped_bases=metagenomic_bases,dataset=EHI_number_api) %>%
  filter(dataset %in% sample_metadata$dataset) %>%
  select(dataset,assembly_length,N50,L50,num_contigs,num_bins,assembly_mapped_bases)
mapping_data <- airtable("tblWDyQmM9rQ9wq57", "appWbHBNLE6iAsMRV") %>% #get base ID from Airtable browser URL
  read_airtable(., fields = c("EHI_sample_static","MAG_mapping_percentage"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
  rename(id=1,dataset=EHI_sample_static) %>%
  filter(dataset %in% sample_metadata$dataset) %>%
  select(dataset,MAG_mapping_percentage)
all_data <- sample_metadata %>%
    mutate(Taxon=factor(Taxon,levels=c("Amphibian","Bird","Mammal","Reptile","Control"))) %>%
    left_join(extraction_data, by=join_by(dataset==dataset)) %>%
    left_join(library_data, by=join_by(dataset==dataset)) %>%
    left_join(indexing_data, by=join_by(dataset==dataset)) %>%
    left_join(preprocessing_data, by=join_by(dataset==dataset)) %>%
    left_join(assembly_data, by=join_by(dataset==dataset)) %>%
    left_join(mapping_data, by=join_by(dataset==dataset))

Laboratory performance differences

DNA yield

Total amount of DNA extracted from the 150 ul subset of the bead-beaten sample.

all_data %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.0f±%.0f", mean(extract), sd(extract))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of total DNA nanograms")
tinytable_hl8pfvwhqvra9ylt6p8q
Mean and standard deviation of total DNA nanograms
Taxon ZYMO DREX EHEX
Amphibian 297±350 340±475 1138±1195
Bird 71±134 11±11 44±45
Mammal 448±432 256±224 867±730
Reptile 102±71 162±134 250±179
Control 0±0 2±3 1±0 
all_data %>%
    ggplot(aes(x=Extraction,y=extract))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="DNA yield (ng)",x="Extraction method")

all_data  %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(extract ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_2uycs2b7jykk4gzw8y86
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 229.412625000 129.0917 1.7771296 18.62337 0.091883533
fixed NA ExtractionDREX -37.448812500 100.5929 -0.3722808 59.99998 0.710995523
fixed NA ExtractionEHEX 345.333000000 100.5929 3.4329752 59.99998 0.001087593
ran_pars Sample sd__(Intercept) 0.003841839 NA NA NA NA
ran_pars Species sd__(Intercept) 373.178642666 NA NA NA NA
ran_pars Residual sd__Observation 348.464099979 NA NA NA NA

Library performance

Number of PCR cycles required for reaching the plateau phase of the indexing PCR.

all_data %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(pcr), sd(pcr))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of optimal number of PCR cycles")
tinytable_ah7lqv2myjiuepjix4ct
Mean and standard deviation of optimal number of PCR cycles
Taxon ZYMO DREX EHEX
Amphibian 14.2±2.1 10.8±3.3 10.2±2.6
Bird 18.0±2.5 18.5±3.9 14.8±2.6
Mammal 11.0±4.0 10.2±1.9 10.3±2.2
Reptile 11.7±4.2 10.3±3.7 12.0±3.5
Control 20.0±2.8 19.8±0.5 22.0±4.2 
all_data %>%
    ggplot(aes(x=Extraction,y=pcr))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="Optimal number of PCR cycles",x="Extraction method")

all_data  %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(pcr ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_ocrk0utygu42jn5reea3
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 13.708333 0.9696194 14.137850 22.50455 1.117002e-12
fixed NA ExtractionDREX -1.250000 0.8896915 -1.404981 60.00000 1.651832e-01
fixed NA ExtractionEHEX -1.875000 0.8896915 -2.107472 60.00000 3.926431e-02
ran_pars Sample sd__(Intercept) 0.000000 NA NA NA NA
ran_pars Species sd__(Intercept) 2.555902 NA NA NA NA
ran_pars Residual sd__Observation 3.081982 NA NA NA NA

Sequencing data

Sequencing depth

#Total GB
all_data %>% summarise(sum(bases_pre_fastp)) / 1000000000
##   sum(bases_pre_fastp)
## 1             429.0952
#Total reads
all_data %>% summarise(sum(bases_pre_fastp)) / 300
##   sum(bases_pre_fastp)
## 1           1430317433
#Average GBs per sample
all_data %>% summarise(mean(bases_pre_fastp)) / 1000000000
##   mean(bases_pre_fastp)
## 1               5.36369
all_data %>% summarise(sd(bases_pre_fastp)) / 1000000000
##   sd(bases_pre_fastp)
## 1            2.703487
#Average reads per sample
all_data %>% summarise(mean(bases_pre_fastp)) / 300
##   mean(bases_pre_fastp)
## 1              17878968
all_data %>% summarise(sd(bases_pre_fastp)) / 300
##   sd(bases_pre_fastp)
## 1             9011624
all_data %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(bases_post_fastp / 1000000000), sd(bases_post_fastp / 1000000000))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of sequencing depth (GB)")
tinytable_9ucewjcfbqptjapdf4rz
Mean and standard deviation of sequencing depth (GB)
Taxon ZYMO DREX EHEX
Amphibian 4.0±1.5 4.2±2.3 4.7±0.4
Bird 3.7±1.9 3.4±2.5 2.7±1.9
Mammal 5.5±3.3 4.6±2.2 3.8±2.5
Reptile 6.1±2.3 5.7±1.4 5.0±1.9
Control 0.0±0.0 0.5±0.6 2.1±2.7 
all_data %>%
    ggplot(aes(x=Extraction,y=bases_pre_fastp))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="DNA yield (ng)",x="Extraction method")

all_data  %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(bases_post_fastp ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_3p91qvdxgmvsdikwikv4
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 4844471692 451697359 10.7250388 3.202405e+01 3.943247e-12
fixed NA ExtractionDREX -386734817 487241483 -0.7937231 3.392332e+20 4.273567e-01
fixed NA ExtractionEHEX -786257695 487241483 -1.6136920 8.048610e+20 1.065942e-01
ran_pars Sample sd__(Intercept) 1195869523 NA NA NA NA
ran_pars Species sd__(Intercept) 555777398 NA NA NA NA
ran_pars Residual sd__Observation 1687854008 NA NA NA NA

Quality-filtering

all_data %>%
    mutate(qf_bases=bases_post_fastp/bases_pre_fastp*100) %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(qf_bases), sd(qf_bases))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of quality-filtered proportion of reads")
tinytable_i318psg6og53jfcz6hkh
Mean and standard deviation of quality-filtered proportion of reads
Taxon ZYMO DREX EHEX
Amphibian 84.6±1.5 91.1±4.8 87.5±4.0
Bird 70.6±16.5 62.5±29.3 66.7±19.3
Mammal 92.2±2.4 89.0±5.1 91.3±1.8
Reptile 89.9±6.6 90.5±7.4 88.3±7.6
Control 3.3±2.3 9.8±11.5 27.5±3.4 
all_data %>%
    mutate(qf_bases=bases_post_fastp/bases_pre_fastp*100) %>%
    ggplot(aes(x=Extraction,y=qf_bases))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="DNA yield (ng)",x="Extraction method")

all_data  %>%
    mutate(qf_bases=bases_post_fastp/bases_pre_fastp*100) %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(qf_bases ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_tl8hlwf4z12ssk6plgkl
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 84.3251792 3.587497 23.5052981 18.58420 2.789261e-15
fixed NA ExtractionDREX -1.0626469 2.798710 -0.3796917 48.00002 7.058490e-01
fixed NA ExtractionEHEX -0.8627821 2.798710 -0.3082785 48.00002 7.592043e-01
ran_pars Sample sd__(Intercept) 6.6948031 NA NA NA NA
ran_pars Species sd__(Intercept) 9.2214285 NA NA NA NA
ran_pars Residual sd__Observation 9.6950144 NA NA NA NA

Genomic data

host_stats <- read_tsv("data/host_stats.tsv")

Read duplicates

host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(duplicates), sd(duplicates))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of fraction of duplicated reads")
tinytable_lqsjzq7t282nxoue55p0
Mean and standard deviation of fraction of duplicated reads
Taxon ZYMO DREX EHEX
Amphibian 0.3±0.2 0.2±0.2 0.2±0.2
Bird 0.8±0.3 0.9±0.1 0.6±0.4
Control 1.0±0.0 0.9±0.0 1.0±0.0
Mammal 0.4±0.4 0.2±0.1 0.2±0.2
Reptile 0.5±0.4 0.3±0.3 0.4±0.4 
host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    ggplot(aes(x=Extraction,y=duplicates))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="Duplication rate",x="Extraction method")

host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(duplicates ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_hskxwrbw58km29kjg4du
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.48796262 0.08012126 6.090302 20.71993 5.105684e-06
fixed NA ExtractionDREX -0.08741146 0.06897524 -1.267287 60.00000 2.099492e-01
fixed NA ExtractionEHEX -0.15230321 0.06897524 -2.208085 60.00000 3.107140e-02
ran_pars Sample sd__(Intercept) 0.00000000 NA NA NA NA
ran_pars Species sd__(Intercept) 0.22019874 NA NA NA NA
ran_pars Residual sd__Observation 0.23893725 NA NA NA NA

Depth of coverage

host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(coverage_depth), sd(coverage_depth))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of fraction of duplicated reads")
tinytable_qeglgch0iopxrn3ejx9k
Mean and standard deviation of fraction of duplicated reads
Taxon ZYMO DREX EHEX
Amphibian 0.0±0.0 0.0±0.0 0.0±0.0
Bird 0.4±0.8 0.6±0.5 0.8±0.8
Control 0.0±0.0 0.0±0.0 0.0±0.0
Mammal 0.7±1.2 0.3±0.4 0.5±1.0
Reptile 0.2±0.3 0.1±0.1 0.1±0.1 
host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    ggplot(aes(x=Extraction,y=coverage_depth))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="Depth of coverage",x="Extraction method")

host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(coverage_depth ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_223fc5mrrmbq9cmm5bsv
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.329125000 0.1331063 2.47264850 20.58494 0.02222872
fixed NA ExtractionDREX -0.089791667 0.1144640 -0.78445305 48.00000 0.43662873
fixed NA ExtractionEHEX -0.002791667 0.1144640 -0.02438903 48.00000 0.98064341
ran_pars Sample sd__(Intercept) 0.408634577 NA NA NA NA
ran_pars Species sd__(Intercept) 0.224731208 NA NA NA NA
ran_pars Residual sd__Observation 0.396515070 NA NA NA NA

Breadth of coverage

host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(coverage_breadth), sd(coverage_breadth))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of depth of host genome coverage")
tinytable_7mzn43y6xrcvqony5mp1
Mean and standard deviation of depth of host genome coverage
Taxon ZYMO DREX EHEX
Amphibian 0.0±0.0 0.0±0.0 0.0±0.0
Bird 0.6±0.5 3.2±4.4 8.9±13.9
Control 0.0±0.0 0.0±0.0 0.0±0.0
Mammal 5.7±5.9 10.2±16.4 15.1±26.4
Reptile 4.9±7.5 3.0±5.4 2.9±3.7 
host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    ggplot(aes(x=Extraction,y=coverage_breadth))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="Breadth of coverage",x="Extraction method")

host_stats %>%
    left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(coverage_breadth ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_ton8kngem6whp7ck6qr9
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 2.799167 2.252826 1.2425133 20.82558 0.22785670
fixed NA ExtractionDREX 1.301250 1.956426 0.6651158 48.00000 0.50915975
fixed NA ExtractionEHEX 3.918750 1.956426 2.0030144 48.00000 0.05084021
ran_pars Sample sd__(Intercept) 7.118462 NA NA NA NA
ran_pars Species sd__(Intercept) 3.549767 NA NA NA NA
ran_pars Residual sd__Observation 6.777259 NA NA NA NA

Metagenomic data

Library complexity

Nonpareil estimate of the metagenomic complexity after removing host DNA.

all_data %>%
    group_by(Taxon,Extraction) %>%
    summarise(value = sprintf("%.1f±%.1f", mean(nonpareil), sd(nonpareil))) %>%
    pivot_wider(names_from = Extraction, values_from = value) %>%
    tt(caption = "Mean and standard deviation of breadth of host genome coverage")
tinytable_spgvkahvr5eq7lc2ld5d
Mean and standard deviation of breadth of host genome coverage
Taxon ZYMO DREX EHEX
Amphibian 0.9±0.1 0.8±0.1 0.8±0.1
Bird 0.9±0.1 1.0±0.0 0.8±0.4
Mammal 0.8±0.2 0.8±0.1 0.7±0.3
Reptile 0.9±0.1 0.9±0.0 0.9±0.1
Control 0.0±0.0 0.5±0.6 0.5±0.7 
all_data %>%
    ggplot(aes(x=Extraction,y=nonpareil))+ 
        geom_boxplot() + 
        facet_grid(. ~ Taxon, scales = "free") +
        labs(y="Nonpareil completeness",x="Extraction method")

all_data  %>%
    filter(Taxon != "Control") %>%
    lmerTest::lmer(nonpareil ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_muhc2owvsupw85vmbwn3
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.885552853 0.03450229 25.6664948 40.89542 7.753071e-27
fixed NA ExtractionDREX 0.005536892 0.04336611 0.1276779 48.00020 8.989373e-01
fixed NA ExtractionEHEX -0.112193866 0.04336611 -2.5871326 48.00020 1.276294e-02
ran_pars Sample sd__(Intercept) 0.059565487 NA NA NA NA
ran_pars Species sd__(Intercept) 0.035030813 NA NA NA NA
ran_pars Residual sd__Observation 0.150224598 NA NA NA NA

Alpha diversity

Variance partitioning metrics are derived from community_analysis.Rmd.

alpha_data <- list.files(path = "results", pattern = "^alpha_.*\\.tsv$", full.names = TRUE) %>%
  map_df(~ read_tsv(.)) %>%
  left_join(sample_metadata,by= join_by(dataset==dataset))
alpha_data %>%
    pivot_longer(!c(dataset,Library,Species,Taxon,Sample,Extraction), names_to = "metric", values_to = "value") %>%
    mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
    ggplot(aes(x=Extraction,y=value))+ 
        geom_boxplot() + 
        facet_grid(metric ~ Taxon, scales = "free")

Richness

alpha_data %>%
    lmerTest::lmer(richness ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_wtyc10ukh0j5e7k6qthi
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 58.43750 12.008400 4.8663854 8.87726 0.0009237868
fixed NA ExtractionDREX 2.18750 4.698926 0.4655319 32.00001 0.6447035294
fixed NA ExtractionEHEX 2.43750 4.698926 0.5187355 32.00001 0.6075141409
ran_pars Sample sd__(Intercept) 15.62413 NA NA NA NA
ran_pars Species sd__(Intercept) 30.71216 NA NA NA NA
ran_pars Residual sd__Observation 13.29057 NA NA NA NA

Neutral

alpha_data %>%
    lmerTest::lmer(neutral ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_v3fqy19chene8hmah3vy
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 29.71528470 5.751141 5.16685008 8.745982 0.0006458923
fixed NA ExtractionDREX -1.59472963 2.085900 -0.76452835 31.999998 0.4501538272
fixed NA ExtractionEHEX -0.04935512 2.085900 -0.02366131 31.999998 0.9812697035
ran_pars Sample sd__(Intercept) 9.00571703 NA NA NA NA
ran_pars Species sd__(Intercept) 14.37531297 NA NA NA NA
ran_pars Residual sd__Observation 5.89981588 NA NA NA NA

Phylogenetic

alpha_data %>%
    lmerTest::lmer(phylogenetic ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
    broom.mixed::tidy() %>%
    tt()
tinytable_g3xk49rgtklm7hymwlx5
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 4.5870074 0.4577012 10.021838 8.589824 4.977993e-06
fixed NA ExtractionDREX 0.2122362 0.1485377 1.428837 32.000064 1.627406e-01
fixed NA ExtractionEHEX 0.2230529 0.1485377 1.501658 32.000064 1.429897e-01
ran_pars Sample sd__(Intercept) 1.1278928 NA NA NA NA
ran_pars Species sd__(Intercept) 0.9754991 NA NA NA NA
ran_pars Residual sd__Observation 0.4201282 NA NA NA NA

Variance partitioning

Variance partitioning metrics are derived from community_analysis.Rmd.

variance_data <- list.files(path = "results", pattern = "^var_.*\\.tsv$", full.names = TRUE) %>%
  map_df(~ read_tsv(.))

variance_data %>% summarise(mean=mean(r2),sd=sd(r2))
## # A tibble: 1 × 2
##    mean     sd
##   <dbl>  <dbl>
## 1 0.102 0.0768
variance_data %>%
    left_join(sample_metadata %>% select(Species,Taxon) %>% unique(),by=join_by(species==Species)) %>%
    mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
    ggplot(aes(x=Taxon,y=r2)) +
    geom_boxplot()+
    ylim(0,1)+
    facet_grid(. ~ metric, scales = "free")

variance_data %>%
    group_by(metric) %>%
    summarise(mean=mean(r2),sd=sd(r2)) %>%
    tt()
tinytable_mj59v9ciiisjtpuentpo
metric mean sd
neutral 0.06201935 0.04744380
phylogenetic 0.07381277 0.06884687
richness 0.16970818 0.06626388